A universal form of slow dynamics in zero-temperature random-field Ising model 
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The zero-temperature Glauber dynamics of the random-field Ising model describes various ubiq- 
uitous phenomena such as avalanches, hysteresis, and related critical phenomena. Here, for a model 
on a random graph with a special initial condition, we derive exactly an evolution equation for an 
order parameter. Through a bifurcation analysis of the obtained equation, we reveal a new class of 
cooperative slow dynamics with the determination of critical exponents. 
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Slow dynamical behaviors caused by cooperative phe- 
nomena are observed in various many-body systems. In 
addition to well-studied examples such as critical slowing 
down [l| , phase ordering kinetics 0] , and slow relaxation 
in glassy systems Q , seemingly different phenomena from 
these examples have also been discovered successively. In 
order to judge whether or not an observed phenomenon 
is qualitatively new, one needs to determine a universal- 
ity class including the phenomenon. In this context, it is 
significant to develop a theoretical method for classifying 
slow dynamics. 

Here, let us recall a standard procedure for classifying 
equilibrium critical phenomena. First, for an order pa- 
rameter m of a mean-field model, a qualitative change in 
the solutions of a self-consistent equation m = T(m) is 
investigated; then, the differences between the results of 
the mean-field model and finite-dimensional systems are 
studied by, for example, a renormalization group method. 
On the basis of this success, an analysis of the dynamics 
of a typical mean-field model is expected to be a first step 
toward determining a universality class of slow dynamics. 

As an example, in the fully connected Ising model with 
Glauber dynamics, an evolution equation for the magne- 
tization, dtm = Q (m), can be derived exactly. The anal- 
ysis of this equation reveals that the critical behavior is 
described by a pitchfork bifurcation in the dynamical sys- 
tem theory 0]. As another example, an evolution equa- 
tion for a time-correlation function and a response func- 
tion was derived exactly for the fully connected spherical 
p-spin glass model 0, Q . The obtained evolution equa- 
tion represents one universality class related to dynamical 
glass transitions. 

The main purpose of this Letter is to present a non- 
trivial class of slow dynamics by exactly deriving an 
evolution equation for an order parameter. The model 
that we consider is the zero-temperature Glauber dy- 
namics of a random-field Ising model, which is a sim- 
ple model for describing various ubiquitous phenomena 
such as avalanches, hysteresis, and related critical phe- 
nomena 0, 0, @, 0, Eli- As a simple, but still non- 
trivial case, we analyze the model on a random graph 
[l2T ] , which is regarded as one type of Bethe lattices 



Thus far, several interesting results on the quasi-static 
properties of the model on Bethe lattices have been ob- 
tained [II, [H, Q3, [HI H, HI. In this Letter, by 



performing the bifurcation analysis of the derived equa- 
tion, we determine the critical exponents characterizing 
singular behaviors of dynamical processes. 

Model: Let G(c, N) be a regular random graph con- 
sisting of N sites, where each site is connected to c sites 
chosen randomly. For a spin variable 0{ = ±1 and a ran- 
dom field h! i on the graph G(c, N), the random- field Ising 
model is defined by the Hamiltonian 
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(1) 



i=l j£Bi i=l 

where Bi represents a set of sites connected to the i site 
and h is a uniform external field. The random field h! i 
obeys a Gaussian distribution DuQi'j) with variance R. 
We collectively denote (ci)i=i and (hi)^ =1 by er and h, 
respectively Let Ui be the number of upward spins in Bi. 
Then, for a given configuration, we express the energy 
increment for the spin flip at i site as — 2eriAi, where 

A i = c-2u i -(h + h' i ). (2) 

The zero-temperature Glauber dynamics is defined as a 
stochastic process in the limit that the temperature tends 
to zero for a standard Glauber dynamics with a finite 
temperature. Specifically, we study a case in which the 
initial condition is given by er^ = —1 for any i. In this 
case, once cr, becomes positive, it never returns. Thus, 
the time evolution rule is expressed by the following sim- 
ple rule: if Cj = — 1 and satisfies A; < 0, the spin flips 
at the rate of 1/tq; otherwise, the transition is forbid- 
den. Note that a. t (t) = — 1 when \(t) > 0, and Aj(t) 
is a non- increasing function of t in each sample [l4j |. In 
the argument below, a probability induced by the prob- 
ability measure for the stochastic time evolution for a 
given realization h is denoted by P h , and the average of 
a quantity X over h is represented by X. 

Order parameter equation: We first note that the lo- 
cal structure of a random graph is the same as a Cay- 
ley tree. In contrast to the case of Cayley trees, a ran- 
dom graph is statistically homogeneous, which simplifies 
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the theoretical analysis. Furthermore, when analyzing 
the model on a random graph in the limit N —> oo, 
we may ignore effects of loops. Even with this assump- 
tion, the theoretical analysis of dynamical behaviors is 
not straightforward, because <jj and <7fc, j,k G Bi, are 
generally correlated. We overcome this difhculty by the 
following three-step approach. 

The first step is to consider a modified system in which 
<7i = —1 is fixed irrespective of the spin configurations. 
We denote a probabilit y in this modifi ed system by Q h . 
We then define q{t) = Q h (<Tj(t) = 1) for j G Bi, where 
q(t) is independent of i and j owing to the statistical 
homogeneity of the random graph. The second step is to 
confirm the fact that any configurations with A j (t) > in 
the original system are realized at time t in the modified 
system as well, provided that the random field and the 
history of a process are identical for the two systems. 
This fact leads to a non-trivial claim that P h (Ai(t) > 0) 
is equal to Q h (Ai(t) > 0). By utilizing this relation, one 
may express P h (Ai(t) > 0) in terms of Q h (<jj(t) = 1). 
The average of this expression over h, with the definition 
p(t) ee P h (A,(t) > 0), leads to 

p® = zZ( c u ) ^n 1 - f 2u h dh ' D R( h ')> 

u=o ^ ' J -°° 

(3) 

where we have employed the statistical independence of 
<7j and <7fc with j, k G Bi in the modified system. The 
expression <j3j) implies that q(t) defined in the modified 
system has a one-to-one correspondence with the quan- 
tity p(t) defined in the original system. The third step 
is to define pit) = Q h (aj = -l,Aj < 0) and r(t) = 
Q h (Aj(t) > 0) for j G Bi. Then, by a procedure similar 
to the derivation of ([3]), we find that dq(t)/dt is equal 
to p(t)/ro. r(t) is also expressed as a function of q(t) 
because r(t) is equal to a probability of Aj(t) > in a 
modified system with <Ji = — 1 and <7j = —1 fixed. Con- 
cretely, we write r{t) = 1 — F(q(t)), where 

F (<?) = E f " ~ 1 ) - ^ U 1^ dh'D R {h'). 

" / Jc-2u-h 

(4) 

By combining a trivial relation p(t) + q(t) + r(t) = 1 with 
these results, we obtain 



u=0 



TO 



dq 
dt 



F(q) - q, 



(5) 



which determines q(t) with the initial condition q(0) = 0. 

Here, for a spin configuration at time t under a 
quenched random field h' t , we define 
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(6) 




FIG. 1: (color online) Phase diagram. Inset: F(q) — q as 
functions of q. R = 1.5. h = 1.0 and h = h c (R). 



limit — > oo. Since is determined by (J3j) with ((5]), 
one may numerically compare p[t) with /5(t) observed in 
the Monte Cairo (MC) simulations. We confirmed that 
these two results coincided with each other within numer- 
ical accuracy. We also numerically checked the validity 
of ([5|) on the basis of the definition of q. As a further 
evidence of the validity of our results, we remark that 
the stationary solution q(oo) satisfies 9(00) = F(q(oo)), 
which is identical to the fixed-point condition of a recur- 
sive equation Z n+ \ = F(Z n ) in the Cayley tree, where 



Z n is a probability at generation n. (See Ref. [14[ for 
the precise definition of Z n .) In the argument below, we 
set To = 1 without loss of generality and we investigate 
the case c = 4 as an example. The obtained results are 
essentially the same for c > 4; however, the behaviors for 
c < 4 are qualitatively different from those for c > 4. 

Bifurcation analysis: We start with the analysis of 
([5]) for R = 1.5 and h = 1.0. The qualitative behavior of 
q(t) is understood from the shape of the graph F(q) — q 
shown in the inset of Fig. [T] There are three zeros, q\, 
q-2, and (73, where < qi < 52 < 93 < 1- Since F(q) > q 
in the interval [0, q\), q(t) — > q\ as t — > 00 with the initial 
condition q(0) = 0. This geometrical structure sustains 
in a region of the parameter space a = (R, h), as shown 
in Fig. [TJ Let qi{a) be the parameter dependence of qi. 
Then, the stable fixed point q\{a) and the unstable fixed 
point qi (a) merge at the solid line, and the stable fixed 
point (73(a) and the unstable fixed point 92(a) merge at 
the dashed line, both of which correspond to saddle-node 
bifurcations The two lines terminate at a critical 
point (R sp , h sp ). Since the trajectories with the initial 
condition <j(0) = do not exhibit a singularity at the 
dashed line, only the bifurcations at the solid line are 
relevant in the present problem. The solid line is called 



21] 



where 6{x) — 1 for x > and 6(x) = otherwise. Due 
to the law of large numbers, p(t) is equal to p(t) in the 



a spinodal line 

Now, we investigate the singular behaviors of slow dy- 
namics near the bifurcation points. We first fix the value 
of R such that < R <C R sp - Then, a saddle- node bifur- 
cation occurs at (R,h c (R)) on the solid line. Let q c (R) 
be the saddle- node point such that q\= qi = q c (R)- We 
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set h = h c (R) + e and q = q c + u. From the graph in the 
inset of Fig. [1] one finds that ([5]) becomes 



du 
~dt 



a Q e + a 2 u 2 + 0(|u| 3 , \eu\), 



(7) 



when | it | <C 1 and |e| <C 1. &o and a 2 are numerical 
constants. Solutions of ([7]) are expressed as a scaling 
form 



u(t) 



,1/2 



u±(\e\ 1/2 t) 



(8) 



when |e| <C 1, where u+ and m_ are e-independent func- 
tions for e > and e < 0, respectively. The result im- 
plies that the characteristic time near q = q c diverges 
as t ~ |e| -1 / 2 . Note that q(t — > oo) = (73 when e > 0. 
Therefore, q(t — > 00) exhibits a discontinuous change, 
and the jump width is given by the distance between 
qi{a) = 92(a) and 53 (a) at e = 0. 

Next, we focus on the dynamical behaviors near the 
critical point (i? sp , h sp ). By substituting q = q c (R sp ) + v 
and (i?, h) = (R sp , h sp ) + (77, e) into ([5]), we obtain 

= c e + cuj« - c 3 w 3 + 0(M 4 , |ew|, |?7w 2 |), (9) 

at 

when |u| <C 1, |e| <C 1, and |^| <C 1. Co, Ci, and C3 are 
numerical constants. System behaviors are classified into 
two types. First, when |?7| 3> |e| 2 ^ 3 , solutions of © are 
expressed as v = \r]\ 1 ^ 2 vi±(\r]\t) when \rj\ <C 1. This scal- 
ing form is identical to that near a pitchfork bifurcation, 
which might be related to conjectures presented in Refs. 
HimH- Second, when \r)\ < |e| 2/3 , which includes 
the case in which R 
expressed as 



R sp is fixed, solutions of 



v=\e\V 5 v 2 ±(\e\h) 



(10) 



when |e| <C 1, with S = 3 and ( = 2/3. The charac- 
teristic time diverges as r ~ |e| near the critical point 
(i? S p, h sp ). In addition to the two scaling forms, one can 
calculate the width of the discontinuous jump of q along 
the spinodal line, which is proportional to (— rf) 1 / 2 near 
the critical point [H, [l6| . 

Finite size fluctuations: In a system with large but 
finite N, fluctuations of p are observed. Their basic char- 
acterization is given by the intensity 



Xp (t) = N((p(t)-(p(t))) 2 )- 



(11) 



where for a quantity X(t) determined by tr(t) and h, 

X(t)) EE Y^P h ^(t)= ( T)X{t)\ c 



{t)=cr- The problem 
here is to determine a singular behavior of Xp under the 
condition that < e « 1 and 77 = 0. Let p c be defined 
by ([3]) with q = q c (R sp ). We then assume 



p(t)-po = A(e,N)F(t/r(e,N)) 



(12) 



near p(t) ~ p c , where A and r are typical values of the 
amplitude and the characteristic time, respectively, and 



F is a time-dependent fluctuating quantity scaled with A 
and r. We further conjecture finite size scaling relations 

A(e,N) = N'^^FAieN 1 /"), 

r(e, N) = N^^FrieN 1 /"), (13) 

where Fa(x) ~ x 1 ' 6 for x ^> 1 and F T (x) ~ x~^ for 
x ^> 1; Fa(x) = const, and F T (x) = const, for x <C 1. 
Here, the exponent v characterizes a cross-over size N c 
between the two regimes as a power-law form iV c ~ e~ v . 
We thus obtain 

Xp (r(e,iV))-7V^F x (e7V 1 /-) ; (14) 

where 7 = v — 2/8. The values of ( and 6 have already 
been determined. We derive the value of v in the follow- 
ing paragraph. 

It is reasonable to assume that the qualitative behav- 
ior of p near the critical point (i? sp , h sp ) is described by 
([9]) with small fluctuation effects due to the finite size 
nature. There are two types of fluctuation effects: one 
from the stochastic time evolution and the other from 
the randomness of h[. The former type is expressed 
by the addition of a weak Gaussian-white noise with a 
noise intensity proportional to 1/N, whereas the latter 
yields a weak quenched disorder of the coefficients of 
©. In particular, e is replaced with e + g/y/~N, where 
j is a time-independent quantity that obeys a Gaus- 
sian distribution. Then, two characteristic sizes are de- 
fined by a balance between the fluctuation effects and 
the deterministic driving force. As the influence of the 
stochastic time evolution rule, the size N s is estimated 
from a dynamical action of the path-integral expression 
/ dt[Nc4(d t v — cq€ + C3V 3 ) 2 + c 5 v 2 }, where the last term 
is a so-called Jacobian term. C4 and C5 are constants. 



the balance among the terms e 2 

4/3 



v 2 /N s 



25l . |26| for a similar 



In fact 

leads to N s ~ e" 4 / 3 . (Sec Refs. 
argument.) On the other hand, the size N q associated 
with the quenched disorder is determined by the balance 

" 2 . Since N s < N„ 



//V q , which leads to 7V q ~ < 
for e«l, the system is dominated by fluctuations when 
N < N q . With this consideration, we conjecture that 
N c = N a ~ e- 2 . That is, v = 2. 



In laboratory and numerical experiments, statistical 
quantities related to the magnetization rh = YliLi a i/N 
may be measured more easily than \p- Since m is not in- 
dependent of p, rh also exhibits singular behaviors near 
the critical point. Concretely, the fluctuation intensity 
x(t) ee N ^(m 2 (i)) — (rh(t)) 2 ^j is characterized by the 
above exponents. In order to confirm this claim, we 
measured x(t) by MC simulations. Then, the charac- 
teristic time is defined as a time t* when x(t) takes a 
maximum value Xm- Our theory predicts i* ~ e -2 / 3 and 
Xm = N 2 / 3 Xm (eN^ 2 ) with xm(x) ~ x~ 4 / 3 for x > 1. 
The numerical results shown in Fig. [5] are consistent 
with the theoretical predictions. 
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FIG. 2: (color online) Time evolution of the magnetization for 
R = 1.781258(~ R sp ) and h = 1.1. Note that h sp = 1. The 
circles are the result obtained by the MC simulations with 
N = 10 5 , and the solid line corresponds to the solution of 
dtm = — (m — 1 + 2p) with ([3]) and ([5]). One may derive this 
equation for m in a straightforward manner. Inset: XmN~ 2 ^ 3 



as functions of eN 1 ^ 2 . The guideline represents 



-4/3 



Concluding remarks: We have derived the exact evo- 
lution equation for the order parameter p describing the 
dynamical behaviors of a random field Ising model on a 
random graph. From this dynamical system, we have de- 
termined the values of the critical exponents: ( = 2/3, 
5 = 3, v = 2, and 7 = 4/3. 

Before ending this Letter, we discuss the critical phe- 
nomena in d-dimcnsional random-field Ising models on 
the basis of our theoretical results. Let v m be the crit- 
ical exponent characterizing the divergence of a corre- 
lation length above an upper-critical dimension d u . By 
assuming a standard diffusion coupling for ^ as an effec- 
tive description of finite-dimensional systems, we expect 
z = C/^m = 2. This leads to v m = 1/3. From an ap- 
plication of a hyper-scaling relation to the upper-critical 
dimension, one also expects the relation v = d n v m [27| . 
This leads to d u = y /y m = 6, which is consistent with 
the previous result [2l|. We then denote the exponents 
characterizing the divergences of time and length scales 
by 1/3 and £3, respectively, for three-dimensional systems. 
These values were estimated a s G — 1.3 and v% ~ 0.7 in 
numerical experiments [22I. [28l.[29|. The theoretical anal- 
ysis of finite dimensional systems will be attempted as 
an extension of the present work; this might be comple- 
mentary to previous studies [U, l3(ij |. Finally, we wish 
to mention that one of the most stimulating studies is 
to discover such a universality class in laboratory experi- 
ments by means of experimental techniques that capture 
spatial structures directly 0, Hi, 31 1 . 
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